% clear all;
% clc;

% ----------------------- Analysis Parameters -----------------------------

pol = -1;
noaxons = 200;

patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

epidur  = 0; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = 0; % 0, -10, -20

patients = 1:5;
nopat = length(patients);

dc_th   = zeros(noaxons,nopat);
dr_th   = zeros(noaxons,nopat);
ra_elec = zeros(1,nopat);

% experimental data
isens_exp_sub = [0.20,0.35,0.70,0.55,0.80];
ipain_exp_sub = [0.50,0.50,1.50,1.40,3.10];
isens_exp_epi = [1.30,19.0,7.00,5.00,9.00];
ipain_exp_epi = [2.20,24.5,19.0,5.00,19.0];    


% ----------------------- Importing Data ----------------------------------

for ii = 1:nopat
    
    patnum = patients(ii);

    % location of electrode
    eletag_epi = 'epi_';
    eletag_sub = 'sub_';
    fdist_elec_epi = cordgeom(patnum,8);
    if(subloc==1)
        fdist_elec_sub = cordgeom(patnum,9);
    elseif(subloc==2)
        fdist_elec_sub = 1-cordgeom(patnum,9);
    else
        return;
    end
        
    loctag_epi = ['p',num2str(patnum),'_d',strrep(num2str(fdist_elec_epi),'.','p'),...
              '_t',strrep(num2str(th_elec),'-','n')];
    loctag_sub = ['p',num2str(patnum),'_d',strrep(num2str(fdist_elec_sub),'.','p'),...
              '_t',strrep(num2str(th_elec),'-','n')];
          
    p_tag = num2str(patnum);

    prefix_epi = ['neg_',eletag_epi,'bp_'];
    suffix_epi = [loctag_epi,'.txt'];
    prefix_sub = ['neg_',eletag_sub,'bp_'];
    suffix_sub = [loctag_sub,'.txt'];
    
    cd(['patient',p_tag,'_thresholds']);

    dcthresh_epi = load([prefix_epi,'dc_act_',suffix_epi]);
    [dcthresh_epi(:,1),dc_si_epi] = sort( dcthresh_epi(:,1)*pol );
    drthresh_epi = load([prefix_epi,'dr_act_',suffix_epi]);
    [drthresh_epi(:,1),dr_si_epi] = sort( drthresh_epi(:,1)*pol );
    
    dcthresh_sub = load([prefix_sub,'dc_act_',suffix_sub]);
    [dcthresh_sub(:,1),dc_si_sub] = sort( dcthresh_sub(:,1)*pol );    
    drthresh_sub = load([prefix_sub,'dr_act_',suffix_sub]);
    [drthresh_sub(:,1),dr_si_sub] = sort( drthresh_sub(:,1)*pol );

    cd ..;

    cd(['patient',num2str(patnum),'_potentials']);
    iunit_epi = load([eletag_epi,'bp_current_',loctag_epi,'.txt']);
    rtiss_epi = 2/iunit_epi/1e3;
    iunit_sub = load([eletag_sub,'bp_current_',loctag_sub,'.txt']);
    rtiss_sub = 2/iunit_sub/1e3;    
    cd ..;
    
    dc_th_epi(:,ii) = dcthresh_epi(:,1);
    dr_th_epi(:,ii) = drthresh_epi(:,1);
    ra_elec_epi(ii) = rtiss_epi;
    dc_th_sub(:,ii) = dcthresh_sub(:,1);
    dr_th_sub(:,ii) = drthresh_sub(:,1);
    ra_elec_sub(ii) = rtiss_sub;  
    
end

RA_epi = ( ra_elec_epi'*ones(1,noaxons) )';
dc_th_epi = dc_th_epi./RA_epi;
dr_th_epi = dr_th_epi./RA_epi;
RA_sub = ( ra_elec_sub'*ones(1,noaxons) )';
dc_th_sub = dc_th_sub./RA_sub;
dr_th_sub = dr_th_sub./RA_sub;

% ----------------------- Processing Data ---------------------------------

isens_model_epi = dc_th_epi(1,:);
ipain_model_epi = dr_th_epi(1,:);
isens_model_sub = dc_th_sub(1,:);
ipain_model_sub = dr_th_sub(1,:);

% ratio_exp   = ipain_exp./isens_exp;
% ratio_model = ipain_model./isens_model;
% ratio_exp
% ratio_model

figure;
hold on;
plot(patients,isens_exp_epi,'ko-','LineWidth',4,'MarkerSize',15);
plot(patients,isens_model_epi,'ko--','LineWidth',4,'MarkerSize',15,...
     'MarkerFaceColor','k');
plot(patients,isens_exp_sub,'ks-','LineWidth',4,'MarkerSize',15);
plot(patients,isens_model_sub,'ks--','LineWidth',4,'MarkerSize',15,...
     'MarkerFaceColor','k');
hold off;
xlabel('Patient Number','FontSize',36);
ylabel('Sensory Threshold (mA)','FontSize',36);
xlim([0.75,5.25]);
ylim([0.1,30]);
set(gca,'FontSize',36,'XTick',1:5,'YScale','Log');
L1 = 'Measured Epidural';
L2 = 'Predicted Epirudal';
L3 = 'Measured Intradural';
L4 = 'Predicted Intradural';
legend(L1,L2,L3,L4,'Location','NE');
